library(ggplot2)
require(data.table)
library("dplyr")
library("ggpubr") 
require(tidyverse) #for easy data manipulation and visualization
require(caret) #for easy machine learning workflow
library(ncdf4)
library(raster)
library("forcats")

wd1 <- "C:\\Users\\yfan\\OneDrive - NORCE\\NFood_replication\\"
case = c("RCP45","RCP45_85LU","RCP85","RCP85_45CO2","RCP85_SAI","RCP85_MSB","RCP85_CCT")
casename = c("RCP45",expression('RCP45'[(SSP5LU)]),"RCP85",expression('RCP85'[(RCP45CO2)]),"RCP85+SAI","RCP85+MSB","RCP85+CCT")


#load data prepared by Data1.GrowingSeasonClimate-Yield.R
load(paste0(wd1,"/data/Rdata/GridSample-SSP2-SSP5-Data-final.Rdata")) 

df.ssp2.all[Year==2090,sum(area),by=.(Scenario,Year)] #this is real yield and area data for RCP45, based on its own area mask (area>1000ha)
df.ssp5.all[Year==2090,sum(area),by=.(Scenario,Year)]
df.ssp5.all[Year==2090,summary(yield),by=.(Scenario,Year)]
df.ssp2.all[Year==2090,summary(yield),by=.(Scenario,Year)] 


#=====draw together fertilizer use and irrigation of SSP2 vs. SSP5
df.lu.samples <- rbind(df.ssp2.all, df.ssp5.all) 
(fertn.crop <- df.lu.samples[,mean(fertn), by=.(Scenario,Crop)]);#(fertn.crop <- df.samples[,max(fertn), by=.(Scenario,Crop)])
df.irrig <- df.lu.samples[Irrig==TRUE, .(yield=weighted.mean(yield, area), area.t=sum(area)), by=.(Year,Scenario,Irrig)] #irrigated means
df.rain <- df.lu.samples[Irrig==FALSE, .(yield=weighted.mean(yield, area), area.t=sum(area)), by=.(Year,Scenario,Irrig)] #irrigated means
df.gl <- df.lu.samples[, .(yield=weighted.mean(yield, area), fertn=weighted.mean(fertn, area), area.t=sum(area)), by=.(Year,Scenario)] #global mean
df.gl <- df.gl[df.irrig, .(AREA=area.t*10^-6, FERTN=fertn, IRRIG=i.area.t/area.t*100 ), by=.EACHI, on=.(Year,Scenario)]
df.gl$Variable=1
df.gl.t <- df.gl[, data.table(t(.SD), keep.rownames=TRUE), .SDcols=-"Variable", by=.(Year,Scenario)]
names(df.gl.t)[3:4] <- c("Variable", "Value") 
df.gl.t[Variable=="FERTN", summary(Value), by=.(Scenario)]
df.gl.t[Variable=="IRRIG", summary(Value), by=.(Scenario)]

#ratio of mean irrigated yield to rainfed yield
df.irrig[Scenario=="SSP2", mean(yield), by=.(Scenario)]$V1/df.rain[Scenario=="SSP2", mean(yield), by=.(Scenario)]$V1


library(facetscales)
library(egg) 
scales_y <- list(
  AREA = scale_y_continuous(limits = c(730, 1020), breaks = seq(750, 950, 100)),
  FERTN = scale_y_continuous(limits = c(3, 18), breaks = seq(5, 15, 5)),
  IRRIG = scale_y_continuous(limits = c(15, 28), breaks = seq(15, 25, 5)))

p <- df.gl.t %>% 
  mutate(Variable = factor(Variable, levels=c("AREA", "FERTN", "IRRIG"))) %>%
  ggplot(aes(Year, Value, group=Scenario)) + geom_line(aes(color = Scenario)) +
  scale_x_continuous(breaks =c(2025, 2050, 2075)) +
  #xlab("") + ylab(expression("   %                gN m"^-2*"            Mha     ")) + 
  xlab("") + ylab(expression(~~~'%'~~~~~~~~~~gN~m^-2~~~~~~~~~~'Mha')) + #use ~ as space, to work with svg format
  #facet_wrap(~Variable, nrow=2, dir="v", scales='free_y', strip.position = "left")
  facet_grid_sc(rows=vars(Variable),  scales = list(y = scales_y)) +
  guides(color = guide_legend( keywidth = 1, keyheight=0.8, label.position = "right", label.hjust = 0, ncol=1))
p <- p + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())+
  #ggtitle(expression('('*bold(e)*') Land use')) + theme(plot.title = element_text(hjust = 0.5,size=15)) +
  ggtitle(expression('Land use')) + theme(plot.title = element_text(hjust = 0.5,size=13)) +
  theme(legend.text = element_text(size=11), legend.title = element_blank(), legend.background = element_blank(), 
        legend.position = c(0.76, 0.77)) + 
  theme(axis.title=element_text(size=11), axis.text = element_text(size=11)) +
  theme(panel.spacing = unit(0, "lines")) +  #reduce facet spacing
  theme(strip.text = element_text(size = 11))#, strip.background = element_blank(), strip.placement = "outside")

tag_facet(p, #tag_facet replace facet strips with in-panel tags
          x = -Inf, y = -Inf, vjust = -6.7, hjust = -0.1,
          open = "", close = "", # style for "(" and ")"
          fontface = 1, size = 4,  family = "sans",# "serif",
          tag_pool = c("Total crop area", "Nitrogen fertilization", "Irrigated proportion"))

ggsave(paste0(wd1,"/figures/main/Fig1e.svg"), device ="svg", units = "in", width = 2.5, height = 3.6, dpi = 600)






#==============supplementary figures for LUC analysis======================
#==============Scatterplot of nitrogen fertilizer and yield
#nitrogen is the major driver of land use change effect
df.ssp2.wm <- df.ssp2.all[, .(yield=weighted.mean(yield, area,na.rm=T), fertn=weighted.mean(fertn, area,na.rm=T), area.t=sum(area,na.rm=T)), by=.(Year,Scenario,Crop)] #
df.ssp2.irrig <- df.ssp2.all[Irrig==TRUE, .(yield=weighted.mean(yield, area,na.rm=T), fertn=weighted.mean(fertn, area,na.rm=T), area.t=sum(area,na.rm=T)), by=.(Year,Scenario,Crop,Irrig)] #
df.ssp2.wm <- df.ssp2.wm[df.ssp2.irrig, f.irrig:=i.area.t/area.t*100, by=.EACHI, on=.(Year,Scenario,Crop)]
df.ssp5.wm <- df.ssp5.all[, .(yield=weighted.mean(yield, area,na.rm=T), fertn=weighted.mean(fertn, area,na.rm=T), area.t=sum(area,na.rm=T)), by=.(Year,Scenario,Crop)]
df.ssp5.irrig <- df.ssp5.all[Irrig==TRUE, .(yield=weighted.mean(yield, area,na.rm=T), fertn=weighted.mean(fertn, area,na.rm=T), area.t=sum(area,na.rm=T)), by=.(Year,Scenario,Crop,Irrig)] #
df.ssp5.wm <- df.ssp5.wm[df.ssp5.irrig, f.irrig:=i.area.t/area.t*100, by=.EACHI, on=.(Year,Scenario,Crop)]

dif.lu.wm <- df.ssp5.wm[df.ssp2.wm, .(Scenario="SSP2-SSP5", yield.dif=log(i.yield)-log(yield), #yield.dif=(i.yield-yield)/yield*100, 
                                      fertn.dif=(i.fertn-fertn)/fertn*100, irrig.dif=(i.f.irrig-f.irrig)/f.irrig*100), by=.EACHI, on=.(Year,Crop)] #


df.ssp5.new <- df.ssp5.all[area>5000 & yield>0.5] 
df.ssp2.new <- df.ssp2.all[area>5000 & yield>0.5] 
summary(df.ssp5.new$yield);summary(df.ssp2.new$yield)
dif.lu.sample <- df.ssp5.new[df.ssp2.new, .(Scenario="SSP2-SSP5", yield.dif=(i.yield-yield)/yield*100, # yield.dif=log(i.yield)-log(yield), #
                                            fertn.dif=(i.fertn-fertn)/fertn*100, area.ssp5=area), by=.EACHI, on=.(Year,Crop,Irrig,Sample)] #
is.na(dif.lu.sample) <- sapply(dif.lu.sample, is.infinite)
(dif.lu.sample[,length(unique(Sample)), by=.(Crop)]);(dif.lu.sample[,summary(yield.dif), by=.(Crop)])
dif.lu.sample <- dif.lu.sample[Crop=="Tropical Corn", Crop:="Corn"]
dif.lu.sample <- dif.lu.sample[Crop=="Tropical Soybean", Crop:="Soybean"]
dif.lu.sample.wm <- dif.lu.sample[,.(yield.dif=weighted.mean(yield.dif, area.ssp5,na.rm=T), fertn.dif=weighted.mean(fertn.dif, area.ssp5,na.rm=T)), by=.(Year,Crop,Irrig)]



palette("default")
library(RColorBrewer)
cbPalette <- brewer.pal(n =8, "Set1")
cbPalette <- c(cbPalette[1:4], "gray2", cbPalette[7])
barplot(height=rep(1,length(cbPalette)),col=cbPalette) #display colors

cropname = c("Maize", "Sugarcane", "Wheat", "Rice", "Soy", "Cotton")
p0 <-  dif.lu.sample.wm[fertn.dif<1000] %>% 
  mutate(Crop = factor(Crop, levels=c("Corn","Sugarcane","Wheat","Rice","Soybean", "Cotton" ), labels=cropname)) %>%
  ggplot(aes(fertn.dif, yield.dif, color = Crop)) + 
  geom_point() + 
  geom_smooth(method=lm) +  # Add linear regression line (by default includes 95% confidence region)
  scale_colour_manual(name = "", values = cbPalette[1:6]) +
  scale_x_continuous("Change in nitrogen fertilizer use (%)")+ #limits =c(0,1000)) +
  scale_y_continuous("Change in yield (%)") #,limits =c(0,1000))
p0 <- p0 + theme_bw() + theme(panel.grid.minor = element_blank()) + #, panel.grid.major = element_blank()
  ggtitle("(a)") + theme(plot.title = element_text(hjust = 0.5), title = element_text(size=9)) +
  theme(strip.text.x = element_text(size = 8)) +
  theme(legend.text = element_text(size=9), legend.title = element_blank(), legend.position = "bottom") + 
  theme(axis.text = element_text(size=8), axis.title=element_text(size=9))

tag_facet(p0, #tag_facet replace facet strips with in-panel tags
          x = -Inf, y=150, #y = -Inf, vjust = -16, hjust = -0.1,
          open = "", close = "", # style for "(" and ")"
          fontface = 1, size = 4,  family = "sans",# "serif",
          tag_pool = c("SSP2 - SSP5"))
ggsave(paste0(wd1,"/figures/supplementary/Fig.S4.svg"), units = "in", width = 6, height = 5, dpi = 600)







